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Abstract. We obtain a simple direct derivation of the differential equation governing 
the entropy flow probability distribution function of a stochastic system first obtained 
by Lebowitz and Spohn. Its solution agrees well with the experimental results of Tietz 
et al [2006 Phys. Rev. Lett. 97 050602]. A trajectory-sampling algorithm allowing 
to evaluate the entropy flow distribution function is introduced and discussed. This 
algorithm turns out to be effective at finite times and in the case of time-dependent 
transition rates, and is successfully applied to an asymmetric simple exclusion process. 
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The concept of entropy is usually associated with probability distributions 
(ensembles) via Gibbs's formula. However, it has been recently emphasized that in 
nonequilibrium systems one can consistently define the entropy production along a single 
trajectory [|1|, |, i, and that the fluctuation theorems |, |, |, |, 0, |, |, |Tg, 0, 0] 
can be interpreted as giving connections between the probability of entropy-generating 
trajectories with respect to that of entropy-annihilating trajectories. Since the entropy 
flow along a given trajectory is experimentally accessible (see, e.g., [|l3l), it is of some 
interest to investigate the distribution function of the entropy production and flow in a 
nonequilibrium system. The distribution of entropy flow in a stochastic system satisfies 
a differential equation derived by Lebowitz and Spohn in ^j. They made use of this 
equation to investigate the large-deviation function of entropy production in the long- 
time limit, thus obtaining a form of the Gallavotti-Cohen symmetry P]. In the present 
paper, we provide a new direct derivation of this differential equation in a general 
nonequilibrium system with discrete states, and show that its solution satisfies more 
general fluctuation relations also at finite times. This derivation highlights its connection 
with the entropy production along a single trajectory. We exploit it to evaluate the 
entropy flow in experiments on an optically driven defect center in diamond [1^. We 
suggest a practical computational scheme to evaluate the distribution of entropy flow. 
We illustrate the feasibility of the method by applying it to the simple asymmetric 
exclusion process (ASEP) [p^ . 

We consider a system with a discrete phase space, whose states are denoted by 
i = 1, . . . ,N. We assume that the evolution of the system is described by a Markovian 
stochastic dynamics 

^ = E mAt)pAt) - w,.m{t)] , (1) 

where pi(t) is the probability that the system is found at state i at time t, and Wij{t) 
is the transition rate from the state j to the state i at time t. We consider a generic 
path uj defined by uj{t) = ik iS tk < t < t^+i, with k = 0,1, . . . , M, with tu+i = t{, and 
define the time-reversed path uj by uj{t) = ik iS t E [tk+i,tk), where t = to + t{ — t. Let 
Q{uj) be defined in terms of the ratio between the probability V{uj) of the forward path 
UJ (conditioned by its initial state io) and the probability V{uj) of the time- reversed path 
UJ, conditioned by its initial state iM = H and subject to the time-reversed protocol 



Qiuj) = - In 



---J2^n "^"^"-^p^ . (2) 

k=l L^ik-iiki^k) \ 

We have assumed that, if Wij{t) > at any time t, one also has Wji{t) > 0. In 
order to interpret this quantity, let us consider the Gibbs entropy of the system 
S{t) = —J2iPiit)lnpi{t) (we take the Boltzmann constant ks = I throughout). 
Using ([^), the time derivative of Ss reads 
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By exploiting the relation Inx < x — 1, one sees that the first sum is non-negative, 
and can thus be interpreted as the entropy production rate 0. The second sum 
defines the entropy 5'f which fiows into the reservoir: dSf/dt = J^i^^j ^ijPj i^ij/^ji)- 
It is easy to verify that, if we define {Q)^ as the average of the quantity Q(a;), 
defined by (^, over all possible paths up to time t, the following equality holds: 
d (Q)^ /dt = dSf/dt, and thus (Q)^^ = ASf = //j dt dSf/dt. Notice that if the detailed 
balance conditions holds for the transition rates Wij{t), and the energy Ei{t) is associated 
to the state i of the system, we have Wji(t)/Wij(t) = exp {[Ei(t) — Ej(t)] /T}, and thus 
Tin [Wji{t) /Wij{t)] represents the heat exchanged with the reservoir in the jump from 
state j to state i. Thus the quantity Q{uj), defined by (@), is the entropy which fiows into 
the reservoir as the system evolves along the path u (see, e.g., 0). We now introduce 
the joint probability distribution (f)i{Q,t), that the system is found at time t in state 
i, having exchanged a total entropy fiow Q with the thermal bath. The entropy which 
fiows into the reservoir as a result of the jump of the system from state j to state i is 
given by Asij = \og[Wji{t) /Wij{t)]. In a small time interval r the variation of (f)i{Q,t) 
is given by 

0i(Q, t + r) ~ 0i(Q, t) + r ^ Wij(j)j{Q - As. 

-Asi 




.n=0 



^{Q- As,,,t)-Wji(Pi{Q,t) 



and thus we obtain the differential equation governing the time evolution of the 
distribution function (j)i{Q,t): 



dUQ,t) 

dt 




.n=0 



-As,,rd^<Pj{Q,t) 



-Wj4i{Q,t) 



By introducing, for each z, the generating function 

iji{\t) = JdQ e^«0,(Q,t), 
and taking into account the expression of Asij, we obtain the master equation 



(5) 



(6) 



dt 



E 



w. 



ijj{\t)-W,ii)i{\t) 



Y,H,,{\)i,,{\t). 



(7) 



This last equation was first derived by Lebowitz and Spohn in 

We now introduce the total probability distribution (j){Q,t) = J2i(pi{Q,t), and the 
total generating function ?/^(A, t) = J2i i'i^, t)- In the case of time-independent transition 
rates Wij, or of transition rates which depend periodically on the time, it can be useful, 
in order to evaluate the distribution function to introduce the large-deviation 

function. In the long-time limit, the generating function ?/'(A,t) is dominated by the 
maximum eigenvalue g{X) of the matrix H(A) = {Hij{X)), which appears in the master 
equation (|^). Therefore, we have, for long times t, 

^{X,t)^exp[tg{X)]. (8) 
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By using the last equation and inverting ( 
the entropy flow in the long time limit as 

dA 



one obtains the probability distribution of 



0(g,t) 



27ri 



(9) 



where A* is the saddle point value implicitly defined by dg/dX\^t = Q/t. If we introduce 
the entropy flow per unit time q = Q/t, we obtain the large-deviation function 



f(q)^g(\*)-X*q 



lim -log0(g,t). 

t^oo t 



(10) 



Note that the functions g{X) and /(g) are Legendre transform of each other, and can be 
then interpreted in terms of path thermodynamics: g{X) can be viewed as a path Gibbs 
free energy, while /(g) is the corresponding Helmholtz free energy. This analogy was 
first pointed out in |]T3|, |TB| where work probability distributions of systems driven out 
of equilibrium were studied. The Gallavotti- Cohen theorem |^ (p{Q) / (p{—Q) = exp(Q) 
is immediately recovered by considering that H(— A + 1) = H"'"(A), where H"'' is the 
transpose of H, and thus has the same eigenvalues as H 

Another remarkable fluctuation relation has been recently proposed by Seifert @]. 
Since the path probability densities V{uj) and V{uj) must be normalized, one obtains 
from (I), for any distribution of initial states p° and any distribution of final states pi 
the following fluctuation relation. 



exp 



Q{uj) + In 




E 
E 

«0>«f 



{tt)=if 

'^(tf)=«0 
Q{to)=i{ 



f 

,QH Pi{_ 
Flo 



1. 



Starting from (|^), this relation is recovered as follows. Let ip'^°{X,t) be the solution 
of (0) with the initial condition ^/'*''(A, t=0) = and let ipj{X,t) = J^igi^fi^^t)- 

The function '?/'j(A=l,t) = 1, Vt, is the solution of (0) with the initial condition 
^^.(A=l,t=0) = 1. Thus, ^ 

= Erfo^?(^=l'^) exp[\n{p,it)/pl 



exp 



JIT]) reads 



«0J 



Y.pj{t)^,{x=i,t) 



(12) 



for any value of t. This result is a generalization of equation (2.20) of ^ and has been 
recently applied to a specific case in [0. 

If the system is characterized by a small number of states, one can explicitly solve 
the equations (|^), and thus obtain the total generating function ipi^Xyt) = J^ii^iiXjt) = 
(e^«) . 

As an example, we consider an optically driven defect center in diamond, which can 
be viewed as a two-state system |T^ . If excited by a red light laser, the defect exhibits 
fluorescence. By superimposing a green light laser, the rate of transition from the 
non-fluorescent to the fluorescent state {W+ in the following) and from the fluorescent 
to the non- fluorescent state (W^), turn out to depend linearly on the green and red 
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light intensity, respectively. In |]T3|], the experimental set-up was such that the rate 
W_ = (21.8 ms)~^ was kept constant, and the rate W+ was modulated according to the 
sinusoidal function W+(t) = Wq [1 + 7 sin(27rt/tm)], with Wq = (15.6ms)"\ 7 = 0.46, 
tm = 50 ms. In the same reference, the histogram of the entropy flux, as defined by (^, 
was measured, over 2000 trajectories of time-length 20 tm- By solving (0) for this system, 
we obtain the generating function ilj{X,t), and using (|^) we obtain the probability 
distribution of the entropy flow (f){Q,t = 20 tm), which is found to agree very well 
with the measured histogram, see figure |l|. 
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Figure 1. Entropy production of a two-state system: a defect center in diamond 
is optically driven from a fluorescent to a non-fluorescent state. The experimental 
histogram is taken from the full line corresponds to the entropy flow distribution 
function (p{Q,t = 20i„i) as obtained from (|^-||). Note that the definition of q in jl^ 
is the negative of ours. 



However this direct approach becomes rapidly impracticable, as the system phase 
space size increases. On the other hand, the evaluation of (j){Q,t) or 4'{X,t) by direct 
simulation of the stochastic process described by the master equation is a highly 
difficult task. Similarly to what happens when one tries to evaluate free energies 
differences by using the Jarzynski equality , dealing with the probability distribution 
of the work of a driven system (see, e.g., |T^, |l6l), the most relevant contributions to 



the average {exp{\Q)) often come from the tails of the distribution of (j){Q,t). 

In ref. fl^, the authors proposed a procedure to evaluate the large deviation 
function g{X) based on biased dynamics and on the parallel evolution of system clones. 
Here we discuss how (^ leads to an alternative scheme for the evaluation of the 
fluctuations in the entropy flow. Differently from what discussed in ref. [0, such a 
scheme allows one to evaluate correctly the function il^{X,t) = (exp(AQ)) at any time, 
and not only in the long time limit, and turns out to be effective also in the case of 
time-dependent transition rates. In the following we adapt to our problem the concept 
of weighted trajectory ensemble, introduced by Sun in and extended by Oberhofer 
et al in where it was successfully applied to the Jarzynski equality. The idea is to 
sample trajectories in the weighted ensemble V{ut) exp [AQ(co't)] rather than successions 
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of single states in the unbiased ensemble. Since ilj{X,t) = J Vujt V{ujt) e-^'-^^^*\ we have 
^ = W).*-(A.4 (13) 

where (. . .);^ is the average in the weighted ensemble V{ut) exp [A(5(ti^t)] /ip{X,t), where 
ijj{X,t) = Z\ is the "partition function" of this ensemble, which will be called the "A- 
ensemble" in the following. The solution of (|1^) thus reads 

(14) 



?/'(A, t) = exp 



dA'(g),, 



While Sun [|I9| proposes a Montecarlo procedure for the sampling of trajectories, so 
as to obtain the average ^Q) ^ the weighted ensemble, in the present paper we propose 
a different procedure which generates trajectories in a suitable entropy-flow weighted 
ensemble. The direct simulation of trajectories in the A-ensemble is hindered by the fact 
that one should already know the exact expression of the ensemble partition function, 
i.e., the function ?/^(A,t), which is the unknown quantity at issue. 

To avoid the problem of the direct evaluation of ^(A, t), following [20], we introduce 
a generic functional of the paths n(u;), and write 

I v^t (g(^t)/n(^,))n(^,)p(^,)e^Q(-') (g/n),,n .... 

^^'>^ jVuJt{l/Yi{uj,))Yi{ut)V{uj,)e^Q^-^) (l/n),_n ' ^ ^ 

where (. . .);^ ^ indicates the average in the new V{uJt)Yi{uJt) exp [XQiut)] ensemble (the 
(A, n) -ensemble in the following). 

In order to simplify the following discussion we consider a stochastic dynamics with 
a discrete small time scale r, such that the jumps between states take place at discrete 
times tk = kr. Within this discrete scheme, we want to write the probability of a path 
in the (A, n)-ensemble. 

The probability of a given path u reads 

where the transition probabilities Kij are defined as Kij = rWij, and Ka = 1 — 
J2j{^{lKji- We now define the new transition probabilities K^j = rWij (Wji/Wij)^, 
and Kii = 1 — Kji, and choose the functional Il{uj), such that 

M 



n(^) = nn.„,,_,(tfc), (17) 



fc=i 



with 



Uij{t) 



I5 if ^ 7^ j ; 

~ (18) 
Knit)/ Knit). iU=J. 

Recalling the definition of Q{uj), (@), we obtain that the probability in the (A, 11)- 
ensemble is given by 

V{uj)Ii{uj) exp [\Q{u)] = %,,,^,K,,_,,,_, ■ ■ ■ K^,,,pI, (19) 

and thus {Q)^ can be evaluated by using (|15|) , for the particular choice of 11, as given by 
(|T8|). Note that (|l^) implies that the one can generate a trajectory in the (A, n)-ensemble 
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by simply simulating the process with the Kij{t) transition probabilities. Furthermore, 
the algorithm here described can be used to evaluate the generating function (0) 
also in the case of time-dependent transition rates. This issue will be addressed in 
a forthcoming paper. 




-2-10123 



A 

Figure 2. Plot of g{X) as obtained by combining (^) and (p^, for the ASEP model. 
The function (?(A) vanishes for A = 0, 1 which corresponds to the normalization 
condition and to Seifert's fluctuation relation (111]), respectively. 
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Figure 3. Histogram of the entropy flow per time unit q, corresponding to 1000 
unbiased trajectories of the ASEP model. Full line: probability distribution function 
(j){q,t{), obtained with the trajectory sampling algorithm. Inset: full line, large- 
deviation function / as a function of the entropy per time unit q, as defined by (|lo|); 
dotted hue: plot of f{q) + q. 



In order to give an example of the application of the approach above described, 
we now consider a nonequilibrium system characterized by a large phase space, and 
evolving according to a stochastic dynamics. Instead of considering systems driven 
out of equihbrium by manipulation as in [jl6|, we study here a steady-state non- 



equilibrium system, namely the asymmetric simple exclusion process (ASEP) |jT^, [22 



it consists in a one- dimensional lattice gas on a lattice of L sites. Each site of the model 
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is either empty or occupied by at most one particle. Each particle can jump into an 
empty nearest neighbor site with transition rates per time unit (rightward) and 
(leftward). The system is kept in a out-of-equilibrium steady state since its first and 
last site are in contact with two particle reservoirs, at densities pa and pb respectively. 
By taking pA > Pb and W+ > W^, one observes a net particle current from the left to 
the right reservoir. We choose the following values of the parameters: L = 100, W+ = 1, 



= 0.75, Pa = 0.75, pb = 0.25, which correspond to the maximum current phase [14 
A direct evaluation of the function g{X) by solving the 2^°° equations (0) is of course 
impossible. We thus apply our trajectory simulation approach to the ASEP model. We 
consider trajectories of time length tf = 5, with the elementary time step r = 0.01: at 
each time the transition probability between two states is given by the transition matrix 
Kij(t). For each value of A we generate Af = 1000 sample trajectories and calculate the 
entropy flow Q, as defined by (0), for each trajectory. Then, for the given value of 
A, by averaging over the Af trajectories, we compute the quantity (Q))^ using ([T5|). 
Finally, by combining @ and (p!4D , we obtain the function g{X) which governs the long 
time behavior of '?/'(A,t). This function is plotted in figure |^. It can be seen that it 
vanishes at A = 0, 1 and is symmetric with respect to A = 1/2. The fact that g{0) = 
corresponds trivially to the normalization condition over all the possible trajectories. 
On the other hand, the fact that the function g vanishes at A = 1, is a non trivial 
result, and corresponds to Seifert's fluctuation theorem (0). The symmetry around 
A = 1/2 corresponds to the Gallavotti-Cohen fluctuation relation. We are now able to 
calculate the large-deviation function /(g) defined in (p!0|), which is plotted in the inset 
of figure ^, together with the expression /(g) +q, which exhibits the symmetry required 
by the Gallavotti-Cohen relation. We check as follows that the quantity /(g) actually 
gives the entropy distribution function 0(g, t) oc exp [t f{q)] for the present model. We 
simulate the unbiased diffusion process (i.e., we use the transition matrix Kij), and 
measure the entropy flow along 1000 trajectories. We then plot the histogram of the 
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measured entropy flow per time unit, together with the function exp [t /(g)], see figure |. 
The agreement between the histogram and the predicted entropy distribution (j){q, t) 
is excellent. Note that it is hopeless to verify the Gallavotti-Cohen and/or the Seifert 
relation (^) by direct inspection of the histogram, since no single point of the histogram 
lies in the range q > 1, where the function f{q) + q has its maximum. This is the same 
problem that one faces when trying to exploit the Jarzynski equality to evaluate the 
free energy of simple microscopic systems . 



Finally, we consider the dynamic states which contribute most to the function g 
for each value of A. We measure the current J(A) of particles that, in the steady state 
of weighted ensembles, jumps to the right (positive current) or to the left (negative 
current), in the unit time. This quantity can be easily measured as biased trajectories 
are generated using the probability Kij. The current J, as measured in both the A and 
(A, n) -ensemble, is plotted in figure as a function of A. It can be clearly seen that, as A 
goes away from the origin (either in the negative or in the positive direction), the current 
differs more and more from its value in the unbiased dynamics (A = 0). Notice that the 
current vanishes at A = 1/2, and that, for larger values of A, J becomes negative, i.e., 
the net motion of the particles is opposite to the density gradient. Thus, the parameter 
A, which is the thermodynamic conjugate of q, selects dynamical trajectories in the same 
way as an external field selects states in an ordinary statistical ensemble. Large values 
of |A| select trajectories which are highly unlikely to appear in an unbiased system. This 
corroborates the interpretation of the large deviation function as a kind of free energy 
in the path thermodynamics [p!5| , |1^ . 

We have thus shown how it is possible to evaluate consistently the probability 
distribution function of the entropy fiow in nonequilibrium systems, by solving the 
differential equation (|^ via the simulation scheme for the (A, n)-ensemble. In this way 
the properties of the trajectories which contribute most to the entropy flow in the regime 
of interest can also be evaluated. 

A.I. is grateful to M. Gianfelice for long and interesting discussions, and to C. 
Mejia-Monasterio for important remarks. 
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